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ABSTRACT 

We present results from a numerical study of the multiphase interstellar medium in sub-Lyman- 
break galaxy protogalactic clumps. Such clumps are abundant at z = 3 and are thought to be a major 
contributor to damped Lya absorption. We model the formation of winds from these clumps and 
show that during star formation episodes they feature outflows with neutral gas velocity widths up to 
several hundred kms^^. Such outflows are consistent with the observed high- velocity dispersion in 
DLAs. In our models thermal energy feedback from winds and supernovae results in efficient outflows 
only when cold ( < SOOi^T), dense ( > lOOMopc^^) clouds are resolved at grid resolution of 12 pc. 
At lower 24 pc resolution the first signs of the multiphase medium are spotted; however, at this low 
resolution thermal injection of feedback energy cannot yet create hot expanding bubbles around star- 
forming regions - instead feedback tends to erase high-density peaks and suppress star formation. At 
12 pc resolution feedback compresses cold clouds, often without disrupting the ongoing star formation; 
at the same time a larger fraction of feedback energy is channeled into low-density bubbles and winds. 
These winds often entrain compact neutral clumps which produce multi-component metal absorption 
lines. 

Subject headings: galaxies: formation — galaxies: kinematics and dynamics — intergalactic medium 



1. INTRODUCTION 

Current numerical galaxy formation models can suc- 
cessfully reproduce some of the properties of damped 
Lya absorbers (DLAs), such as the lower end (-/Vhi < 

• lO^i-S cm^^) of the column density distribution and the 
total incidence rate (Pontzen et al. 2008; Razoumov et al. 
2007), the distribution of metals, and the slope of the 
relation between metallicity and low-ion velocity width 
which appears to originate in the mass-metallicity re- 
lation in the models (Pontzen et al. 2008). On the 

\ other hand, simulations tend to overpredict the num- 
ber of DLAs with A'hi ^ 10^^'^ cm~^ and systematically 
produce fewer high-velocity systems. Most such systems 
feature multiple components in their absorption line pro- 
files, but unfortunately one cannot map these compo- 
nents from the velocity space to real space to identify the 
absorption regions and constrain the mechanism produc- 
ing such high velocities. 

In general, the velocity dispersion of neutral gas clouds 
can come either from the gravitational infall in the pro- 
cess of hierarchical buildup of galaxies, in the form of ran- 
dom velocities of protogalactic clumps (Haehnelt et al. 
1998), or from feedback from stellar winds and super- 
novae (SNe) (Schaye 2001). In fairly massive 10^^ 
halos at z = 3 as much as 20 — 30% of gas by mass can 
be in the cold phase surviving the infall (Razoumov et al. 
2007). The corresponding Wdrc ^ 250kms~^ can account 
for part of the observed neutral gas velocity dispersion. 
However, more massive halos are rare at z = 3, and the 
fraction of cold gas drops sharply in > 10^^ M© halos, 
leaving us in search of other mechanisms to produce high 
velocities. 

Galactic winds driven by the feedback energy from 
stellar winds and SNe are an obvious candidate (Schaye 
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2001). Star-forming Lyman-break galaxies (LBGs) at 
2 ~ 3 show evidence for large-scale outflows with typical 
velocities of hundreds kms~^ (Pettini et al. 1998, 2001). 
In fact, with a simple semi-analytical model McDonald 
& Miralda-Escude (1999) showed that feedback at the 
rate 1.8 x 10^° erg yr^^ per 10^^ M© of halo dark matter 
mass added to the velocity dispersion of neutral clouds 
inside virialized halos works out perfectly to explain the 
observed DLA kinematics. However, this energy trans- 
fer takes place on pc scales currently inaccessible to cos- 
mological models. Moreover, the inability of numerical 
galaxy formation models to capture physics on such small 
scales has led to a number of predicaments, the most fa- 
mous of which is the overcooling problem, accompanied 
by the excessive loss of angular momentum in simulated 
galactic disks. 

This classical problem (Katz 1992) has been somewhat 
alleviated in recent years (Thacker & Couchman 2000; 
Sommer-Larsen et al. 2003) as it was realized that feed- 
back from young stars can be very efficient at keeping 
gas in a diluted state preventing it from rapid collapse 
and conversion into stars. However, even galaxy mod- 
els at a sub-kpc (0.1 — Ikpc) resolution cannot capture 
propagation of supernova blastwaves into the interstellar 
medium (ISM), as the injected thermal energy is radi- 
ated away very quickly before it can be converted into 
kinetic energy. The reason is very simple: the mass of a 
resolution element to which the feedback energy is sup- 
plied is usually several orders of magnitude larger than 
the typical mass of a SN ejecta. Therefore, the tempera- 
ture and expansion velocity of the post-shock regions are 
greatly underestimated, and so is the cooling time which 
scales as oc 

above 10^ K (Dalla Vecchia & Schaye 

2008). 

Cosmological simulations must then turn to ad-hoc as- 
sumptions about the role of stellar feedback at scales 
below their resolution limit. Two types of solutions have 
been popular. The first one is suppressing radiative 
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cooling in the feedback regions for the duration of the 
starburst (Mori et al. 1997; Thacker & Couchman 2000; 
Sommer-Larsen et al. 2003; Stinson et al. 2006), to allow 
more efficient conversion of feedback energy into hydro- 
dynamical expansion. This approach leads to more re- 
alistic simulated galaxies that correctly reproduce many 
of the observed properties of present-day disk galaxies 
and have only with a small deficiency in angular mo- 
mentum. On the other hand due to lack of resolution 
the feedback energy is supplied to a fairly large mass, 
and the outflow velocities are usually underestimated. 
This method is known to produce "puffy" galaxies that 
cannot reproduce the high-end tail of the observed DLA 
velocity dispersion (Razoumov et al. 2007). Since these 
galaxies extend to larger radii, their outer regions may 
pick part of the velocity dispersion of the local galaxy 
group, so that they have slightly less severe kinemat- 
ics problem (Pontzcn et al. 2008) than similar-resolution 
models which do not suppress radiative cooling. 

The second widespread approach is to use kinetic feed- 
back instead of thermal feedback (Navarro & White 
1993; Springcl & Hcrnquist 2003; Dalla Vccchia & Schayc 
2008), usually implemented in particle-based simula- 
tions. Although there are several variations of this 
method, the basic idea is to give a velocity kick to a 
small fraction of gas particles near the star-forming re- 
gions, adjusting the mass loading and velocities to re- 
produce observations. Some of the recipes turn off hy- 
drodynamical interaction of the wind particles with the 
gas to obtain large-scale outflows, while other stress the 
importance of such interaction to create hot bubbles in 
the ISM and develop galactic fountains (Dalla Vecchia & 
Schaye 2008). 

A popular method to circumvent some of the resolu- 
tion problems that can be combined with either of the 
above two approaches is to use a sub-resolution multi- 
phase model that describes analytically growth of cold 
clouds embedded in a hot intercloud medium, star for- 
mation (SF) in these clouds, feedback and cloud evapo- 
ration (Yepes et al. 1997; Springel & Hernquist 2003). In 
such models SF and feedback are self-regulating. How- 
ever, different phases are not dynamically separated from 
each other, and therefore by itself such model cannot re- 
sult in outflows. 

At the other end of the resolution spectrum, detailed 
models of small patches of galactic disks, usually in the 
context of the Milky Way galaxy, provide sufficient reso- 
lution to study turbulent ISM stirred by SN explosions. 
Such models resolve hot bubbles driven by individual 
SNe, fragmentation of shells created by these bubbles, 
and the structure of cold dense clouds on pc scales (e.g., 
Joung & Mac Low 2006). High SF rates in such mod- 
els naturally lead to galactic outflows, galactic fountains 
rising to several kpc away from the midplane, and shell 
fragments raining back onto the disk as intermediate- 
velocity cold clouds (Joung & Mac Low 2006). Such 
high-resolution simulations can in principle be used to 
develop subgrid models of stellar energy feedback for a 
given SF rate in cosmological simulations, although to 
the best of our knowledge currently there are no subgrid 
models in the literature that separate the dynamics of 
different components of the ISM. 

In the past few years it has become possible to extend 
such high-resolution 3D models to entire galactic disks. 



albeit at a lower spatial resolution. Tasker & Bryan 
(2006) used adaptive mesh reflnement (AMR) models 
to study the multiphase ISM in a quiescent Milky Way- 
sized disk galaxy. They employ two SF prescriptions, one 
based on cosmological simulations, with low SF threshold 
(0.02 cm~^) and low efficiency, and the other one with a 
much higher threshold (lO'' cm^'') and a high efficiency. 
Their highest numerical resolution is 25 pc which is a typ- 
ical size of giant molecular clouds; they include cooling 
to 300 K and later add photoelectric heating (Tasker & 
Bryan 2008). Their models reproduce a multiphase ISM 
with most of the mass in cold, dense clouds, while SN 
feedback drives gas out of the plane of the galaxy, but 
most of it eventually falls back on the disk. All of their 
models reproduce the slope of the observed relation be- 
tween the SF rate and the gas surface density, on both 
global and local scales, although the high-density thresh- 
old models tend to produce more intermittent outflows 
and occasionally triggered SF in the outer disk. 

Saitoh et al. (2008) carried out SPH simulations of 
an isolated gas disk with 10^ — lO'' particles to study 
the effect of various SF prescriptions on the structure of 
the ISM. Similar to Tasker & Bryan (2006, 2008), they 
test both a cosmological (0.1 cm~'^) and a high-density 
(100 cm~^) SF thresholds, but also vary the SF efficiency 
esF- Only the high-density threshold models could repro- 
duce the complex multiphase structure of the gas disk, 
regardless of the value of esF- In these runs the SF rates 
depend on the modeled global gas flow from intermediate 
densities to the actual sites of SF rather then the actual 
prescribed SF efficiency. On the other hand, the low- 
density threshold models produce thicker and smoother 
disks with SF rates highly sensitive to the chosen value 
of esF- Therefore, they conclude, the use of a high SF 
threshold will avoid uncertainties in the SF models. 

Ceverino fc Klypin (2007) developed a realistic pre- 
scription for modeling feedback formulating conditions 
under which simulations would resolve the formation of 
hot bubbles in the multiphase ISM. Such bubble can 
only be created if simulations resolve cold dense clouds in 
which feedback occurs - only in these clouds heating can 
exceed cooling so that a large fraction of feedback energy 
is converted into hydrodynamical expansion, creating the 
hot phase and driving winds. Ceverino & Klypin (2007) 
also added heating by massive binary systems which are 
ejected from molecular clouds when one of the compo- 
nents becomes a SN. These "runaway stars" carry energy 
very efficiently away from high-density regions, eventu- 
ally exploding as SNe in lower density environments and 
thus facilitating the formation of the hot gas component 
even in low-resolution cosmological models. 

In this paper the approach of Ceverino & Klypin (2007) 
is used to study the formation of winds in high-redshift 
protogalactic clumps responsible for damped Lya ab- 
sorption. We show that a brief episode of SF in a sub- 
LBG galaxy that creates a multiphase medium can also 
drive winds with neutral gas velocity dispersions up to 
several hundred km s~^ . If a substantial fraction of ^ = 3 
protogalactic clumps form such winds at any given time, 
these winds can explain the observed DLA kinematics 
(Schaye 2001). Our models resolve the effect of massive 
stars in protogalactic clumps. Although the spatial reso- 
lution of this study (12 pc) is not sufficient to follow the 
details of shell fragmentation or turbulent interactions, 
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Fig. 1. — HI column density map of the most massive halo in 
model Ml in Razoumov et al. (2007) at 2 = 3 (see text for details). 
The projection is approximately 100 kpc on a side. 



it is argued that this resolution is adequate for model- 
ing the multiphase medium for the purposes of comput- 
ing galaxy formation and launching galactic winds in the 
cosmological context. 

2. MODELS 

2.1. Peak cross-section of DLA absorption 

The mass range of halos that are the main contributors 
to the total DLA line density is still debated. Pontzen 
et al. (2008) argue that the main contribution comes 
from halos in the mass range 10^ — 10^"'^ M©. Lower- 
mass halos have much smaller absorption cross-sections 
due to heating by the ultraviolet background (UVB), 
while the number of halos with masses above 10^^ Mq 
drops sharply. Pontzen et al. (2008) point out that their 
peak at ~ 10^° Mq is probably related to their particu- 
lar feedback implementation in which cooling is turned 
off to reproduce the blastwave solution. On the other 
hand, Nagamine et al. (2007) see a peak of DLA ab- 
sorption shift to higher halo masses with the increased 
wind strength, reaching 10^^ Mq in the "strong wind" 
model. 

In our earlier cosmological DLA models (Razoumov 
et al. 2007), the most common DLA absorbers are ha- 
los in the range 10^" — 10^^'^ Mq. At the lower end of 
this range, absorption is typically dominated by a sin- 
gle galaxy in the halo, while in halos with masses above 
~ 10^^ Mq DLAs are commonly associated with one of 
several protogalactic clumps with the average gas clump 
mass of ~ (1 — 2) X 10^ Mq (Fig. 1). Our goal is to model 
formation of winds in one of these clumps as it undergoes 
a starburst. For the sake of simplicity, we adopt a fixed 
gravitational potential with Mhaio = 3 x IQ-'^^ Mq and a 
low /disk = Afdisk/Afhaio = 0.005, but wc arguc that our 
results are equally applicable to lower-mass (10^° Mq) 
halos with /disk ~ 0.05 - 0.10%. 



2.2. Initial setup and grids 

All simulations in this paper were performed using the 
AMR hydrodynamical code ENZO (O'Shea et al. 2004). 
The computational domain is a 3D periodic box 100 kpc 
on a side covered with a 64^ root grid and up to seven 
levels of refinement corresponding to 12 pc spatial resolu- 
tion. A fixed spherical Navarro-Frenk- White DM profile 
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is assumed at the center of the volume, where x = cr/rvir, 
the concentration parameter c = 12, and Mhaio = 
3 X 10^^ Mq. The initial distribution of gas follows an 
isothermal disk with a temperature of 10'' K and a den- 
sity 



p(r,z) = poe-'-/'-"sech2 (^^^ 



(2) 



(Tasker & Bryan 2006). For all models we take zq — 50 pc 
and ro = 800 pc to approximate a compact star-forming 
disk at z = 3. The total mass of the gas disk is then 



Mdisk = /diskMiaio = S-n-parlzo- 



(3) 



Since our disks are at least 30 times smaller than the 
side of the box, using a small 100 kpc volume is a rea- 
sonable approximation; in addition, small perturbations 
from nearby clumps should be always expected. To re- 
solve the initial disk configuration, inside the central 
(20kpc)'^ region a hierarchy of six centered nested grids 
is set, with the maximum initial resolution of 24 pc. We 
start with the uniform temperature T — 10* K, resolving 
the Jeans length everywhere in the disk. During evolu- 
tion we refine adaptively by the local Jeans length re- 
quiring that it must be resolved by at least 16 cells at 
all times, which is four times better than the Truelove 
criterion (Truelove et al. 1997). 

A cooling function in the temperature range T — 
300 — 10^ K is assumed, with heating by the ionizing 
UVB from Razoumov et al. (2006) with self-shielding 
above 0.01 cm~'^. Since DLAs at z = 3 have fairly high 
metallicities at times exceeding solar, for simplicity solar 
metallicity is assumed throughout all calculations. 

2.3. Star formation and feedback 

Star formation is modeled with discrete stellar parti- 
cles that represent a population of stars born in the same 
cell roughly at the same time and assumed to have the 
same velocity vector in later evolution. In two of the 
three runs presented in this paper, we adopt the mini- 
mum stellar particle mass M^^aiin = 100 Mq; our stellar 
particles can have any mass above M^ min if a sufficient 
amount of gas in a given cell satisfies the following SF 
criteria. A stellar particle is always created at the finest 
local level of refinement, in cells in which (1) the total 
gas density exceeds the threshold psF, and (2) the mass 
of the gas is larger than the local Jeans mass. In other 
words, stars will be formed only in cells in which the lo- 
cal Jeans length is unresolved which with our refinement 
criterion is possible only at the highest AMR level. If 
a cell is marked as a candidate for SF, we compute the 
mass of stars it would form with the given efficiency esF 
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over the local dynamical time tdyn, and scale that mass 
to the local timestep At. If /3(A.T)'^esFAt/tdyn exceeds 
M*,min, a stellar particle is created, and the correspond- 
ing mass is removed from the gas component. Since our 
minimum stellar particle mass is very low and would al- 
low us to record individual core-collapse SN events, we 
adopt instantaneous conversion of gas into stars with esF, 
unlike, e.g., in Tasker & Bryan (2008) where the actual 
SF and feedback associated with each stellar particle are 
continuous over the dynamical timescale. 

Over its lifetime every stellar particle injects feedback 
energy into the thermal energy of the gas. We use a 
prescription similar to that of Ceverino & Klypin (2007) 
to include feedback by both stellar winds and type II 
SNe. Stellar winds supply energy at a constant rate of 
3.88 X lO'^-^ ergs^i for nh = 40 Myrs after creation 
of the stellar particle which corresponds to conversion of 
77en = 2.72 X 10~^ of the rest-mass energy of newly cre- 
ated stars into feedback energy. In addition, during the 
last 10% of the 40 Myr feedback phase, SNe contribute 
10^^ erg per every 100 M© of the initial stellar particle 
mass. This energy is added to the thermal energy of 
the cell hosting each stellar particle, in discrete 10^^ erg 
events spread uniformly over the 4 Myrs interval. In our 
starburst model most stellar particle masses fall into the 
range 100 — 200 Mq, therefore our simulations begin to 
resolve environmental effect of individual massive stars 
exploding as SNe. 

In addition to energy, winds and SNe also return mass 
and metals to the ISM. In our models, energy and mass 
release into the ISM is strictly synchronized to avoid 
putting too much energy into regions which have been 
cleared by winds and/or earlier SNe. We assume that 
Vmass = 0-25 of the total mass that goes into stars is 
ejected back into the ISM via winds and SNe. It can be 
shown that the maximum sound speed in hot bubbles in 
simulations in which mass and energy are simultaneously 
released into the same cell cannot exceed 

/ \ 1/2 

Cs ^ ^ c« 1800kms-\ (4) 

\ ^mass / 

where c is the speed of light, corresponding to the maxi- 
mum temperature of ~ 4 x 10® K. Actual temperatures in 
hot bubbles are somewhat lower in the range 10^ — 10® K, 
largely due to expansion and the work performed to com- 
press the ambient medium, and to a much lesser degree 
due to cooling in the bubble itself. With outflow ve- 
locities added to Cg, the Courant-Friedrichs-Lewy (CFL) 
condition sets the shortest timesteps in our models to 
several thousand years. 

In all our runs the SF efficiency esF 0.3 is assumed. 
Many authors have found that the exact value of esF 
has little impact on the mean SF rates (Stinson et al. 
2006, e.g.), as long as it is in the range from 0.05 to 1 
(see their Fig. 14). Moreover, the true efficiency of SF, 
i.e. the fraction of gas that is eventually converted into 
stars in dense clouds should be determined by an inter- 
play of various processes starting from the hydrodynam- 
ical timescale of gas supply to the star-forming regions. 
Saitoh et al. (2008) have found that at sufficiently high 
PSF the SF rates are effectively set by the timescale of 
cold gas supply from reservoirs (nn = Icm^'') to the 
star-forming regions {nu ^ 100 cm~^) which in their cal- 



TABLE 1 

Simulation parameters: (1) highest level of refinement 

'max, (2) SF THRESHOLD psF, (3) MINIMUM STELLAR PARTICLE 
MASS M^^inin- 



model 


^max 


PSf/(M0Pc-3) 


M*,niin/M0 


Al 


7 


158 


100 


A2 


7 


158 


1000 


A3 


6 


25 


100 



culations is about five times longer than the local dy- 
namical timescale in the star-forming regions. 
In this paper we are using a set of three simulations 

listed in Table 1: a high- resolution starburst model Al, 
a high- resolution quiescent disk model A2, and a low- 
resolution model A3 for which psF was adjusted to pro- 
duce the highest SF rate. The high-resolution models 
used 7 levels of AMR in the (10 kpc)^ x 6 kpc region cen- 
tered on the disk resulting in 12 pc grid resolution. The 
low-resolution model employed 6 levels of refinement in 
the same region corresponding to 24 pc spatial resolution. 

We ran the quiescent model to estimate the effect 
of star formation on the structure of the ISM and on 
galactic wind kinematics. There are two ways to re- 
duce the SF rate in our models: increase the SF den- 
sity threshold pgp or increase the minimum stellar par- 
ticle mass M»_i„in. Note that in our setup M^^min is not 
independent of psF- A star particle is formed only if 
p{Ax)^esFAt/tdyn in a cell exceeds M^^y^in- For the fidu- 
cial value PSF = 158M0pc~'^ and esF = 0.3, the gas 
mass in a cell allowed to form stars is 

^-■- = ^^»^^( i58m'pc-3 )(20^)- 

Setting M,^min = 100 Mq would result in immediate SF 
once the density exceeds psF, whereas using a much 
higher value would delay SF until more gas accumu- 
lates in the cell. For the quiescent disk model we use 
M^*,min = 1000 M0. This prescription ultimately results 
in conversion of approximately the same amount of gas 
into stars, but over a several times longer period, and 
produces a very different ISM morphology and much 
weaker winds. 

3. RESULTS 

3.1. Global ISM morphology 

Without cooling our disks would be marginally 
Toomre-unstable. Adding cooling leads to rapid gas ac- 
cumulation near the galactic midplane and its subsequent 
fragmentation into cold clumps and warm interclump 
material. For the initial central disk density 10~^^gcm3 
and temperature 10* K the cooling time is of order of few 
years leading to gas collapse of the disk onto the mid- 
plane on the timescale of ~ 20 Myrs. Soon thereafter 
first cold clumps form in which SF begins. By 50 Myrs 
high-resolution models start developing a complex ISM 
morphology characterized by dense clouds and filaments 
separated by warm (10, 000 — 20, 000 K) medium seen in 
many simulations (e.g., Wada & Norman 2001), and to 
a lesser degree in the low-resolution model. In model Al 
ample gas supply quickly leads to a starburst starting 
at t « 50 Myrs and lasting ~ 80 Myrs (Fig. 2) in which 
~ 20% of the gas in the disk is converted into stars. 
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In the quiescent disk model A2 tliere is no single star- 
burst phase; the first significant episode of gas conversion 
into stars takes place well past t = 100 Myrs, with inter- 
mittent SF throughout the entire run. It is interesting 
that by the end of simulation A2 at t « 420 Myrs ap- 
proximately the same total stellar mass (3 x 10* Mq) is 
accumulated, although its effect on the underlying gas 
distribution will be completely different. 

The morphology of the ISM is clearly affected by feed- 
back from SF as can be seen from the surface density 
maps in Fig. 3 at t = 119 Myrs corresponding to the 
end of the starburst phase in model Al. By this time 
in the starburst model a much larger amount of mass 
and energy have been injected through feedback into the 
lower-density gas. The result is a much larger role of 
pressure confinement of cold clouds in model Al, as op- 
posed to more gravitational confinement in the quiescent 
model. The mass distribution appears to be smoother 
in the starburst model, with a lower density contrast be- 
tween the clumps and the voids (Fig. 4). Also evident in 
model Al is a more pronounced gas accumulation near 
the galactic center, and a violent stripping of the outer 
gas regions of the disk by feedback waves. The latter 
process depends, of course, on the gas mass of the outer 
disk which in turn is determined by the cosmic accretion 
which we do not compute in our current models. 

In our high-resolution models dense clouds are contin- 
uously being formed and destroyed by self-gravity, dif- 
ferential rotation, feedback from SF inside the clouds, 
and interaction with feedback waves coming from nearby 
star-forming regions. Any single cloud usually survives 
only for a fraction of its galactic orbital revolution, in 
other words, from few Myrs to few tens of Myrs. This 
is consistent with many estimates of the giant molecu- 
lar cloud (GMC) lifetimes in the Milky Way galaxy, al- 
though we do not resolve the scales and processes taking 
place inside these clouds. 

3.2. Conditions in sinT,ulated SF regions 

We are interested in modeling conditions in the star- 
forming regions that facilitate launching of galactic winds 
from thermal feedback only, without suppression of cool- 
ing. We will here review a set of criteria necessary to 
model winds and expanding hot bubbles in the ISM. First 
and foremost, heating by a SN must lead to a sharp rise 
in the gas temperature that would drive the hot bubble 
expansion without rapid cooling. In other words, during 
all stages of bubble expansion heating must dominate 
over radiative cooling. For the early stages, this con- 
dition was elegantly formulated in Ceverino & Klypin 
(2007) (see their eq. 5); using our SF threshold, we will 
write it as 



r ^ 7.8 X lO^^ergs-^M" 



A 



10 



22 erg cm'^ s 



158M0PC-3 
1 



(6) 



(7) 



where /o» is the spatial density of young stars, and p*/p 
is expected to be in the range 0.1 — 1. If the gas tem- 
perature is aroimd 10^ K, cooling A ~ 10~22 erg cm"^ s^^, 
and heating from SNe cannot counterbalance cooling in 
any moderate overdensity. On the other hand, at very 
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Fig. 2. — Total SF rates in all three runs sampled at 50 kyr time 
intervals. 



low temperatures (~ 100 K) cooling is much less efficient 

(A ~ 10~25 ergcm^ s"-'^), and even at high star-forming 
cloud densities feedback may be able to heat the gas. 
Therefore, Ceverino & Klypin (2007) argue, it is crucial 
to include cooling to ~ 100 K to resolve the cold phase in 
order to heat up the gas via SN feedback. Once hydro- 
dynamical expansion of the feedback region begins, gas 
flows out, the mass ratio p^/p increases, assisting further 
heating and expansion. 

When a SN injects energy into the ISM, the resulting 
pressure in the hot bubble greatly exceeds the surround- 
ing pressure. Provided that the energy is not quickly ra- 
diated away, the bubble expands only if it is not confined 
by self-gravity. This second condition was formulated in 
Ceverino & Klypin (2007) in terms of the pressure dif- 
ference between the bubble and the surrounding gas 



. ^ 4ttG , .o 



(8) 



where r is the radius of the bubble, and p is the ambient 
gas density. Using the ideal gas equation of state and our 
SF density threshold, we can rewrite eq. 8 to obtain the 
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Fig. 3. — Disk surface density of models Al, A2, A3 at t 
119 Myrs. 
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Fig. 4. — PDF of gas density in the disk, defined as the fraction 
of the volume per unit logarithm density interval (thick solid line), 
in all three models at t = 119 Myrs. Only cells in the r < lOkpc, 
\z\ < 100 pc region are shown. All three components of the ISM — 
cold (T < lO^K, dotted line), warm (10^ K < T < 10* K, dash- 
dotted line), and hot (T > 10* K, dashed line) - are clearly visible. 
Vertical dashed lines show the SF thresholds. 

minimum resolution necessary to model the expansion of 
the HII regions against self-gravity 
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Since our models start to resolve individual SNe, the 
mass of a resolution element should be small enough in 
order for it to get heated by the typical ~ 10^^ erg ex- 
plosion energy. A single SN explosion in a dense cloud 
may have a hydrodynamical impact only if it can heat 
its host cell to high (10® — 10^ K) temperatures of an ex- 
panding hot bubble. In other words, the energy input 
per SN should then exceed 



(10) 



2/ir7iH 

which gives us an estimate of the minimum resolution 
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necessary to heat up the host ceh 



Aa; ^ 2.49 pc 



EsN 



158Mopc- 



1/3 



-1/3 



T 



106 K 



-1/3 



(11) 



(12) 



At first glance, this constraint requires much higher res- 
ohition than the self-gravity condition (cq. 9). Fortu- 
nately, type II SN explosions have a 30 — 40 Myr delay 
after the initial starburst, and many of tliciii explode in 
environments which have been previously cleared by stel- 
lar winds and neighbouring SNe. Even more importantly, 
the lifetimes of individual cold clouds arc usually in the 
range from few Myrs to few tens of Myrs. By the time a 
stellar particle hosts a SN explosion, its cloud of origin is 
very likely to have been destroyed, and the SN energy is 
released into a p <C psF environment. In addition, stel- 
lar particles may have non- negligible intrinsic velocities - 
traveling even a small 5kms~^ velocity for 35 Myrs will 
take a particle 180 pc away from its birthplace. There- 
fore, the ISM densities in which type II SN explosions 
take place are likely to be several orders of magnitude 
smaller than psF making the constraint in eq. 12 much 
less demanding. 

Once conversion of feedback energy into hydrodynami- 
cal expansion becomes efficient, lack of spatial resolution 
can present an additional problem. If the density con- 
trast in the ISM is not resolved in the simulation, heating 
and hydrodynamical stirring of the star-forming clouds 
might erase high-density peaks bringing the ongoing SF 
to a halt. In other words, SF/feedback can be too self- 
regulating at low resolution. Since we do not know a pri- 
ori the amount of clumping in the ISM at 2: = 3, perhaps 
the most reliable way to reduce this eff'ect is to compare 
the SF rates at various resolutions. 

3.3. Neutral gas kinematics in quasar absorption lines 

Fig. 5 shows the p — T diagram of the entire simulation 
volume in all three runs at t = 119 Myrs. In cells host- 
ing stellar particles, stellar winds and SNe return both 
mass and energy. Assuming that winds remove all ambi- 
ent gas, the minimum density in such cells can be easily 
estimated from the mass loss rate of each star particle of 
mass M« during its feedback stage 



Plo 



0.25M, Ax 1 
40 Myrs 2ufiow Aa;3 

^ / ^flow \ -1 
VlOOOkms-V 



8.6 X 10"^ cm"^ 



(13) 



100 Up 



Ax 
12 pc 



where Wflow is the fiducial outflow speed, and Ax is the 
cell size. These feedback regions can be easily seen in the 
starburst model in Fig. 5 at T ~ 10^ K, and to a much 
lesser degree in the quiescent model. 

Since we model isolated systems without external ac- 
cretion, winds from the disk do not experience any ram 
pressure of the infalling material and can be stopped only 
by gravity and collision with the gas previously blown 
off the disk. The wind speeds often exceed several hun- 
dred kms~^, consequently a large fraction of the vol- 
ume is quickly flUed by a low-density (~ 10"*^ cm~^) gas 




-2 
log r„/cm~^ 

Fig. 5. — P — T diagram of the entire simulation volume in each 
run at 119 Myrs. Each plot is divided into 200'^ colls, and each 
cell is colored by its mass fraction (the darker means a higher mass 
fraction) . 

(Fig. 5). Since we model a large (100 kpc) computational 
volume, the mass fraction locked in this low-density com- 
ponent is significant, especially in the quiescent model 
A2. In the starburst model Al the density of the wind is 
clearly much higher (Fig. 5). Can such dense winds from 
starburst environments account for the wide absorption 
line profiles seen in DLAs? 

To answer this question, we constructed a set of low- 
ionization metal line spectra. At each time output, we 
projected 200 random lines of sight within 3 kpc of the 
center of each disk and calculated absorption line pro- 
files of an unsaturated low-ion transition along sight-lines 
with HI column density above 10^°'^ cm~^. For each such 
profile we calculated a line width wgo corresponding to 
90% of the total optical depth of all components in the 
line. Note that this diagnostic measures the neutral gas 
velocity dispersion, not the typical outfiow velocities, and 
is dominated by clouds with large optical depths its more 
detailed discussion and the comparison to the equivalent 
width can be found in Prochaska et al. (2008). Fig. 6 
shows the distribution of ugo widths for each model as a 
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Fig. 6. — vgo absorption velocity widths of low-ionization lines 
in three models as a function of time. The solid line in each panel 
shows the median velocity. The dashed horizontal line indicates 
the observed 90kms~^ DLA median velocity. 

function of time, along with the median value. Although 
these velocities cannot be compared to DLAs statistics 
directly, since we do not have a cosmological sample and 
do not account for gas accretion which would regulate 
SF episodes, wc note that the typical absorption veloci- 
ties are much higher in the star burst model. None of the 
cosmological DLA models can reproduce the observed 
median {vgo) w 90kms~^, with the actual value of (ugo) 
of order 40 — SOkms"-'^ in Razoumov et al. (2007) and 
close to 60kms^^ in (Pontzcn et al. 2008). Note that in 
cosmological models the simulated widths are also sensi- 
tive to the DLA cross-sections, as, for instance, the ve- 
locity dispersion in "puffy" galaxies with extended radial 
profiles will have an additional weight. 

We argue that if a substantial fraction of clumps along 
the quasar line of sight in a host halo experience an active 
starburst, it could result in a much larger velocity dis- 
persion possibly explaining the observed incidence rate 
of high-velocity DLAs. The fraction of such active star- 
forming galaxies is poorly constrained and can be only 
computed from the cosmic gas infall rate onto individual 



systems in cosmological simulations. 

Here instead we focus on individual systems. In the 
starburst model {vqq) approximately matches the ob- 
served value for 60 Myrs, whereas the quiescent model 
has {vgo) ~ 50kms~^ for most of the disk evolution. The 
low-resolution model features a delayed and much weaker 
SF resulting in a brief episode during which (wgo) reaches 
the observed value. Fig. 7 shows the time evolution of the 
volume fractions of cold, warm and hot gas in the plane of 
the disk. After very rapid cooling the early {t < 50 Myrs) 
evolution is marked by separation of gas into warm and 
cold components. Shortly thereafter feedback gives rise 
to the hot component, and subsequent evolution of the 
disk is characterized by recurrent episodes of gas heating 
and cooling. The "depth" of these episodes is higher in 
the low-resolution model, in which relatively low-density 
cold chimps are more susceptible to feedback. At high 
resolution the relative change of the volume fraction of 
cold and hot gas is visibly reduced, as the number of cold 
star-forming clumps increases, and so does the density in 
individual clumps. The higher SF rate results in a larger 
volume filling fraction of the hot gas, at the same time 
driving more energetic winds from the disk. 

We can see the formation of hot galactic chimneys driv- 
ing the outflows in the vertical slice in Fig. 8. The vi- 
sual structure of the outflows is very different from the 
single-source models of winds from high-redshift dwarf 
galaxies (e.g., Fujita et al. 2004), more resembling the 
high-resolution multiphase models of Ceverino & Klypin 
(2007); Wada (2008). Our current spatial resolution is 
not yet sufScient to model instabilities in the shells cre- 
ated by hot bubbles (Ferrara & Ricotti 2006) or even 
follow these shells farther away from the disk. 

Fig. 9 shows the range of velocities and densities in the 
wind in the starburst model. Ionized, low-density wind 
moves with velocities up to several thousand kms~^, 
whereas neutral gas exhibits velocities of few hundred 
kms~^. Cold gas absorption comes from \z\ < 3kpc; 
there is a clear asymmetry above and below the disk, ex- 
plained by the fact that feedback starts in a fairly small 
number of cold clumps. This asymmetry can be also seen 
in the HI column density maps in Fig. 10. It is important 
to remember that since our models do not account for in- 
teraction of winds and shells with the infalling material, 
the extent of the HI absorbing regions might change in 
more realistic models with cosmic infall. 

4. CONCLUSIONS 

We use high-resolution hydrodynamic; simulations of 
isolated z = 3 protogalactic clumps to show for the first 
time that the high-end tail of the DLA neutral gas veloc- 
ity width distribution can be naturally produced if a sub- 
stantial fraction of clumps in 10^° — few x 10^^ Mq halos 
experience a starburst event at the rate ^ 5 Mq yr^""^ for 
a sustained period of few tens of Myrs. Such starbursts 
would produce a multiphase ISM in which the volume 
fraction of the hot gas exceeds 50% in the midplane of 
the clump. Hot bubbles created by feedback expand both 
in the horizontal direction to fill the intcrclump material, 
and vertically to form chimney-like structures which lead 
to high- velocity winds. Shells and neutral gas fragments 
embedded in such winds 1 — 3kpc above the disk give 
rise to neutral gas absorption with the median velocity 
{vqo) ~ 90kms~^ lasting roughly for the entire duration 
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Fig. 7.— Volume fraction of cold (T < 10^ K, blue), warm 
(10^ - 3 X 10"' K, green), and hot (T > 3 X 10"* K, red) gas in the 
galactic midplane (\z\ < 50 pc) for all three models. 



of tlie starburst. 

Similar to other multiphase disk models, the cold com- 
ponent of the ISM is distributed in a complex network of 
filamentary structures confined by hot bubbles and voids. 
Inside these filaments dense clouds form as gravitational 
instabilities the growth of which is assisted by the ex- 
ternal pressure. Once their growth begins, the clouds 
proceed to form stars very quickly. Only a small frac- 
tion of each cloud (~ 1 — 2% by mass) is converted into 
stars, as the clouds are being continuously destroyed by 
stellar winds from inside, interaction with shells pushed 
by rapidly expanding nearby hot bubbles, and collisions 
with other clouds. In our simulation clouds are transient 
objects constantly exchanging mass and energy with the 
cold shells and filaments, as well as with the hot gas, 
and rarely surviving as discrete entities for longer than 
a fraction of the rotational period. 

In this paper we chose the typical clump masses and 
initial conditions representative of protogalactic environ- 
ments at z = 3. We expect similar winds to arise in 
grid-based cosmological simulations at ^ 10 pc spatial 
resolution. Although none of the current cosmological 
models have such resolution, we argue that using AMR 
techniques to zoom in only on those protogalactic clumps 
that have a high gas infall rate will allow to obtain galac- 
tic winds from thermal feedback in the cosmological con- 
text, without suppressing cooling or using kinetic feed- 
back. 
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Fig. 9. — Outflow kinematics in model Al at t = 119 Myrs. 
Top panel: scatter plot of gas density vs. the vertical velocity 
component. Lower panel: height above/below the disk vs. the 
vertical velocity component for all njj > 0.01 cm~^ cells. 
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Fig. 10. — HI column density of model Al at t = 119 Myrs. 



